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I . INTRODUCTION 


1 2 

In previous works * the authors have described a sequential estimator 
(called the J-Adaptive estimator) which tracks the state of a dynamical system 
with a simplified system model by simultaneously tracking a forcing term which 
approximates the unmodeled system dynamics. This estimator has been very 
successfully simulated in satellite orbit determination problems , where the 
unmodeled system dynamics result from unmodeled earth oblateness accelerations. 
The J-Adaptive estimator requires, however, the specification of a statistical 
model for the unmodeled system dynamics. 

The present study extends the J-Adaptive Estimator to include the simultan- 
eous estimation of the statistics of the unmodeled system accelerations, thus 
completely automating the estimation process. Simulations in satellite orbit 
determination demonstrate the effectiveness of the resulting estimator. 

The report is organized as follows. Section II reviews the J-Adaptive 
estimation concepts previously reported in detail. The algorithm for estimating 
the unmodeled acceleration noise statistics is presented in Section III, which 
also gives the full estimator equations. Section IV presents the simulation 
results in detail. Finally, Section VII presents the conclusions of this study 
and recommendations for further work. 
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11 . THE J-ADAFTIVE ESTIMATOR 

The basic concepts behind the J-Adaptive estimator have been reported 
1 2 

previously. ’ • In essence, the real dynamical system 

R = f(R,R), (1) 

where R is the satellite position vector (R = V) , is modeled in the estima- 
tor as ’ 

R - f m (R,R) + S(t)u, (2) 

so that S(t)u approximates 

f(R,R) - f m (R,R), (3) 

which is the modeling error made in the estimator dynamics. S(t) is a 
specified "symmetry" matrix (not contained in the work reported in Ref [1]), 
and u is the unmodeled acceleration vector. 

The J-Adaptive estimator sequentially tracks the state R,V and the 
unmodeled acceleration vector u from tracking (or other) data. In the 
estimator, the unmodeled acceleration u is modeled statistically as a random 
polynomial in time (see Ref [1] and Section III). The equations for the basic 
J-Adaptive estimator are given in Eqs (24-29) of Ref [1]. The inclusion of the 
symmetry matrix S modifies the prediction equation for the state estimate x 
(first of Eqs (24) of Ref [1}), and the definitions of the ty, ^ and 
matrices of Eqs (14) of Ref [13, which now become 
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III. ESTIMATION OF NOISE STATISTICS 


The statistical model taken for the random (unmodeled) acceleration vector 
u ‘in Ref [1] is of the form 

u(k+l) = u(k) + u(k)x, x = t fc+1 - t k> (5) 

where 

u(k+l) = u(k) (6) 

and the covariance matrix of u is specified a ? priori to be some (diagonal) 

matrix U. . which is not allowed to decrease as a result of the estimation 
uu 

process (as In a "consider" filter mode). On the other hand, u and u are 
estimated (unlike the "consider" filter mode). This is the basic J-Adaptive 
filter. 

Specification of U in the above filter requires some engineering judge- 
ment and experimentation in a given problem, and has been found to be informa- 
tion rate dependent. It is clearly desirable to automate the process of selection 
of ' U by making it data dependent, or adaptive. Such adaptive estimation of 
U (or some other statistic in a random model for the acceleration u) is the 
subject of the present report. 

The basic approach taken to the estimation of the statistics of the 
acceleration u is the adaptive filtering approach of Refs [3,4]. In essence, 
at each measurement time k, that value of the statistic is selected which 
produces the most likely average measurement residual. This adaptive process 
for the selection of. the statistics of u will be made explicit below within 
the context of the J-Adaptive filter. 

The statistical model ultimately selected for u is given by 

u(k+I) = u(k) + Tw(k) (7) 

where {w(k) } is a 3-vector, zero-mean, white Gaussian sequence with identically 
distributed components and covariance matrix 



( 8 ) 


e{w(k)w T (£) } = q(k)I 6^ , 


where q(k) is to be estimated from the measurements. The model for u given 

in Eqs (5-6) was abandoned in favor of the model given in Eqn (7) for the 

following reasons. The filter with the model in Eqs (5-6) contains correlation 

matrices (correlation between the state and u) and (correlation 

between u and u) . These correlation matrices force the state covariance matrix 

P and the u covariance matrix U . C. and U . are in turn forced by 

uu u uu 

U. . . When U. . fluctuates randomly, as it does when it is estimated from the 
uu uu 

measurements, it introduces numerical inaccuracies in the correlations and 

U . . These in turn cause ill- conditioning of P and U , sometimes to the 
uu uu 

point of loss of positive definiteness. 

The complete dynamical model selected is therefore given by 


x(k+l) = <f> (x(k), Su(t», u(t) = u(k) + (t-t.)w(k) 
m k 

u(k+l) = u(k) + xw(k) 


(9) 


where x is the state vector (R,V) and <f> m is the mapping defined by Eqn (2). 
The measurement model is given by 


y i (k) - h i (x(k)) + v i (k), i = 1,N 


( 10 ) 


where y^ are the scalar measurements available at time k, and (v^(k)} 
are independent, scalar, zero-mean,' white Gaussian noise sequences with 


e{v i (k)vt(f)> = R,,(k)6, 


^Li^kA * 


(ID 


The measurement model is written in vector* form as 


y(k) - h(x(k)) + v(k) , 

where 

y T = E^i — y^l * h T = [hi ••• _hj, R = diag(R ii ). 


( 12 ) 


(13) 
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Now define the sensitivity matrices 


$(k+l»k) 63 


r 3x(k+iji = 

r 36 -r 

T m 

|_3k(k) J 

L^(k) J 


¥(k+l,k) 



r(k+i,k) * 


3x(k+l j~ 
3w(k) 


A 



M(k) - 


M x (k) 

v k > 


.9ji(xCk)_l 
_ 3x(k) 


( 14 ) 


and the covariance and correlation matrices 

P(k| *.) = e{[]c(k) - x(k|«3 g(k) - i(k|*0 T } 

c u (kU)' - E { [i(k) - i(kU)3 E(k) - u(k|<on T } as) 

U(k|i) = e{[u(k) - u(k|j03 [u(k) - u(k|4)] T } 

where P(k|'Jl) is the state estimation error covariance matrix (covariance of 
errors in x(k|jt)) at time t^, given all measurement up to and including 

time t 5 . C and U have similar meaning. Then the usual extended Kalman 

Xr u 

filter^ with variance q(k) known and specified is given below. 


Prediction (in time) is given by 

x(k+l|k) = 4> m (x(k|k), Su(t|k)) 
u(k+l |k) - u(k|k) 


(16) 
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P'(k+l|k) *P(k|k)* T + *C (k|k)Y T + yc^(k|k)* T + ¥U(k|k)^ T 

P (k+1 | k) = P"(k+l|k) + q(k)rr T 
C u (k+l|k) * $C u (k|k) + YU(k[k) + q(k)xr 
U(k+l|k) = U(k| k) + q(k)T 2 X 

where the argument (k+l,k) has been omitted from the matrices $>, 
The (measurement) update is given by 

x(k|k) = x(k|k-l) + K (k) £(k) - h(i(k|k-l)3 

A 

u(k|k) = u(k|k-l) + K^(k) [y(k) - h(x(k |k-l))3 

P(k|k) = P(k|k-1) - K (k)M(k)P(k|k-l) 

C u (k|k) = C u (k|k-1) - K x (k)M(k)C u (k|k-l) 

D(k|k) - D(k|k-1) - K u (k)M(k)C u (k|k-I) 

where 

K (k) = P(kik-l)M T (k)Y _1 (k) 

K u (k) = J(k|k-l)M T (k)Y' 1 (k) 

and where 

Y(k) » W(k)P(k|k-l)M T (k) + R(k) 


( 17 ) 


y and T. 


(18) 


(19) 


( 20 ) 


( 21 ) 
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Now che variance q(k) will be estimated via the techniques of Refs (3, 4] 
To that end, define the average normalized predicted residual at time k+1, 


'k+1 N 


! N ~ \(*(k+l|k)Ll 

M ' L 


Z=1 


4<k + n 


( 22 ) 


It is easy to compute 


e( r ^ +1 > = e(r^ +1 |q(k)=0} + dq(k) 


(23) 


where 


e{r^ +1 |q( k )=°} 


, N N 

n + “T l l 

2=1 m=l 


M 5 (k+1 )P" (k+1 | k)M^(k+l) 


4 (k+i) c< k+i) 


(24) 


. N N M 0 (k+l)rrV(k+l) 

l ^ ^ ■ m 

N 2 2= 1 m=l R^(k-Hl) R^(k+1) 


(25) 


3,4 


Then the most probable q(k) based on r^ + ^ is given by 

■ 

i ['k+i - “ e{ ? k+il q(k)E0} J’ if p° sl£ive 


q(k) = 


, otherwise 


with a = 1. 


(26) 


The estimate q(k) of q(k) utilized is an exponentially age-weighted 
average of q(k) generated from 


c fc = ac^j + 1, - 0, a < 1, 

q(k) = ^ q(k-l) + -f- q(k), q(-D - 0 . 
C k C k • 


(27) 
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The q(k) from Eqs (27) is used for q(k) in the filter equations (17). 

To complete the estimator equations, the symmetry matrix S in Eqs (14) 
and (16) must be specified. A matrix S used with success is given by 





IV. DETAILS OF SIMULATIONS 


The simulations presented in the next section feature a circular, polar 
orbit of approximately 1000 km altitude. The "real" system force model consists 
,'of the two-body model, plus the earth oblateness model given in Appendix A, 
plus a mascon (mass concentration) buried at a depth of 100 km. The mascon 
acceleration is represented as 




R - R 


m 


R - R 


m 1 


+ 



(29) 


where k is of course the satellite position vector. R is the mascon 

32 ‘ m _8 

position vector, and p = 0.02 km /sec tso that p - 5x10 p, where p 

m m 

is the gravitational parameter of the earth] . In summary, the real system is 
represented by 


R = - + AR + AR (30) 

r 3 o m 

To complete the specification of the system constants (see Appendix A) , the 
rotation rate of the earth and the eccentricity of the earth are, respectively, 

<o “ 7. 2921159x10" 5 rad/sec, e =0.081813336 . (31) 

Tracking (range and range rate) is simulated from several tracking stations 
to give continuous , and often overlapping, station coverage, sampled every 5 sec. 
Tracking begins at the north pole and usually continues for 50 minutes, which 
is slightly less than one-half the orbital period. The error models for range 
and range rate measurements are simplified, consisting only of additive white 
Gaussian noise. Measurement noise standard deviations vary, and are specified 
for each simulation. 

Estimator dynamical models are variously two-body, and two-body plus 
approximate and oblateness accelerations. Precise estimator models 

used are specified in each simulation. 
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V. SIMULATION RESULTS • 

The first set of simulations to be described feature the J-Adaptive esti- 
mator with estimated noise variance q (J-A-q estimator) with a two-body 

A 

dynamical model. The J-A-q estimator is also compared with the J-Adaptive 
estimator with engineered a priori statistic (J-A-U estimator), also 

with a two-body dynamical model. Note that, with a two-body model, closed form 
prediction equations are available, and no numerical integration is required. 

• Figures 1-7 show the performance of the J-A-q estimator with a two-body 

dynamical model, when the measurement noise standard deviations for range and 

range rate are o^ = 3m, a. = 1 cm/sec, respectively. As indicated in Section 

IV, continuous and often overlapping tracking station coverage is available. 

Figures 1-3 show the tracking of the unmodeled acceleration. The solid lines 

represent the actual unmodeled acceleration (oblateness and mascon) , while the 

dots are the estimates Su. The z-component of the unmodeled acceleration 

-5 -4 2 

dominates, being of the order of 10 (initally 10 ) km/sec , and is there- 

fore tracked most precisely. The x and y acceleration component estimates 
are more noisy, but are clearly tracked. Note that there is an initial transient 
lasting about 300 sec or so before the J-A-q estimator locks onto the unmodeled 
acceleration. Figures 4 and 5 show the normalized measurement residuals (nor- 
malized by their noise standard deviations) for range and range rate. The 
initial transient described above exhibits itself in the range rate residuals, 
producing initially, rather large residuals. This transient dies out, however, 
and the measurement residuals then appear quite white and are within the 
measurement noise level. The absence of residuals after about 2000 sec of 
tracking is merely a programming problem in the computer-generated plots; the 
residuals continue to exhibit appropriate statistical properties for the balance 
of 1000 sec. Figures 6 and 7 show the actual state estimation errors in 
position and velocity; namely | R— R | and j V— V | . It is seen that after an 
initial transient the state estimation errors are quite comparable with the 
measurement noise. That -is, position is estimated to about 3m and velocity to 
about 1-5 cra/sec. 

Let us discuss the initial - 300 sec transient which is observed in 
Figures 1-7. This transient is, in part, due to the transient in u ; that is, 

BLANK NOT FILMED 
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the state cannot be estimated precisely until the unmodeled acceleration is. 

This is not the whole story, however, Note that the range rate normalized 
residuals, while sometimes as large as 5a initially, are random and cannot 
explain the large initial state estimation errors, particularly in velocity. 

It is noted that the drop in state errors at about 300 sec coincides with 
acquisition of the satellite by a second tracking station. It is also noted 
that the position and velocity estimation errors are largely in the x and y 
components; the z component errors are quire small. This is simply due to the 
station-satellite geometry. Thus from the first station, x, y, x and y are 
not as observable as z and z are. This geometric observability clearly has 
an impact on the estimation of Su as well; and also on the estimation of q. 
q guarantees residuals consistent with a projection of the state covariance 
matrix on the available measurement space; not on the unobservable subspace of 
state space. This last observation was well evidenced in a simulation with 
range rate measurements only. In that case, range rate residuals were just as 
nice as in the present simulation, while state errors were not quite as good 
during the first 300 sec of tracking, q compensated only for the state sub- 
space it could see. This leads to the observation that while additional meas- 
urement types (e.g. angles) may not enhance a standard estimation process 
(e.g. batch), they may well enhance the J-A-q estimator. 

The precision of tracking of the uiunodeled acceleration can be improved 
with more precise measurements. This can be seen in Figure 8, which shows the 
x acceleration component estimate generated by the J-A-q estimator with 

o. = 0.3m and o. - 0.1 cm/sec. Compare this with Figure 1. 

P P . 

The J-A-U estimator was simulated on the case described above (a - 3m, 

_P 

a. - 1 cm/sec). This was easy to do once the J-A-q run was made; U was set 

P A 

at the steady-state value of ql. Figure 9 shows the tracking of the x 
acceleration component. Comparing this with Figure 1, improvement in tracking 
is seen, although the J-A-U estimates definitely have deterministic oscillations 
and are not as desirable from this point of view. Figure 10 shows that the 
J-A-U estimator produces somewhat nicer measurement residuals than the J-A-q 
estimator (compare with Figure 5). But a stand-alone J-A-U run would require 
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engineering the U statistics which would normally take several simulations* 

In contrast to this, all J-A-q estimator runs are completely automatic, and are 

untouched by human hand . 


The next set of simulations to be described feature the J-A-q estimator 
with an improved dynamical model. In addition to the two-body term, the esti- 
mator contains the J 2 and terms with approximate values of and J^; 

namely, = 1 

simulations is to determine how well the estimator might identify the higher 
order geopotential. Figure 11 shows the tracking of the z acceleration com- 
ponent with Op - 3m, = 1 cm/sec. It is observed that the lower accelera- 

tion levels cannot be tracked with such noisy measurements. However, Figures 
12 and 13 show that the state is still estimated to the data precision; the 
undetected accelerations are inconsequential relative to the overall noise 
levels. Figure 14 shows the tracking of the z acceleration component with 

a = 0.3m, a. = 0.1 cm/sec. • Tracking is improved (compare with Figure 11), 

P P 

but the extremely low unmodeled accelerations are still undetected. Things 
of course look much better on a linear scale (see Figure 15) . Tracking can be 
further improved with - 3 cm, a. » 0.01 cm/sec, as can be seen in Figure 16. 
When better instruments are developed we will show additional simulations. 


.083 x 10 


-3 


4 


J 3 = -2.55 


x 10 


■ 6 4 


The objective of these 


The next simulations feature an unmodeled mascon. That is, the only model 
error is due to the mascon. As seen in Figure 17, the J-A-q estimator (unaided 
by human hand) cannot track this extremely small unmodeled mascon acceleration. 

As can be seen in Figure 18, however, a carefully engineered J-A-U estimator 
does detect the mascon, albeit with some lag and overshoot which eventually 
damps out. To accomplish this, however, we need = 1 cm and a ^ = 0.01 cm/sec. 


Limited simulations were performed on a single station pass. This is the 

first station of the earlier simulations (recall the earlier station-satellite 

geometry discussion) . The measurement noise standard deviations for range and 

range rate are = 3m , ■ a ^ - 1 cm/sec, respectively. Table 1 summarizes 

the position and velocity errors at the end of the pass for several cases. 

4 

First a standard extended Kalman filter run was made, where the filter model 
was perfect (full oblateness). The lack of good observability in the x and y 
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coordinate directions is clearly seen in Table 1 (although it should be noted 

4 

that the filter was not iterated, either locally or globally; some small 
improvement might be expected with iteration.) Next a totally uncompensated, 
standard extended Kalman filter run was made, where the filter model contained 
the two-body term and the J ^ and terms with the approximate values of 

and used above. It is seen in Table 1 that this uncompensated filter 

is diverging, as expected. A J-A-q estimator, with the two-body and the same 
approximate ^ and J^ terms was. completely unsuccessful (started to diverge). 
This is due to the initial transient, discussed earlier, where the u estimates 
are poor; the J-A-q estimator is not yet locked onto the unmodeled accelera- 
tion. However, improved results are obtained (shown in the last line of Table 1) , 
when the estimated u is not used for prediction. This might be termed an 
"adaptive consider mode" of the estimator. As seen in Table 1, these results 
(third line) are comparable to the perfect model run (first line) in the x 
and y coordinate directions. z-coordinate errors are larger, however; the 
J^A-q adaptive consider mode filter is putting too much uncertainty in the 
z-direction. 

These single station pass simulations indicate that when observability is 
marginal, automatic adaptive estimation is of questionable value. Great care 
must be exercise in adding uncertainty to a selected subspace of the state space; 
so as not to destroy the information content of the data. 
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Table 1 


SINGLE STATION PASS ERRORS 



Position Errors (m) 

Velocity Errors (m/sec) 

A 

x-x 

A 

y-y 

A 

z~z 

• 

• * 

x-x 


A 

z-z 

Perfect Model 
Extended Kalman 

146. 

-171. 

- 9.8 

0.22 

■ 

0.003 

Uncompensated 
Extended Kalman 

558. 

-646. 

-68.6 

0.80 

-0.95 

-0.070 


150. 

-118.. 

-335. . 

H 

-0.07 

-1.34 






























VI. NEW TECHNOLOGY 


The work under this contract consisted of the development of recursive 
estimation techniques for application to satellite orbit determination. The 
previously developed J-Adaptive estimator was modified and extended to include 
the geometry of random (unmodeled) accelerations, and the adaptive estimation 
of noise statistics in a stochastic model for these accelerations. The 
algorithms developed were tested by computer simulation on satellite orbit 
problems with unmodeled earth oblateness and mascon accelerations of various 
types. 

The algorithm developed as a result of this effort is applicable to satel- 
lite orbit determination when excellent, continuous tracking coverage is 
available, as may be the case in satellite- to-satellite tracking. 

Frequent reviews and a final survey for new technology were performed. 
While the extensions to the J-Adaptive estimator developed represent new 
mathematical techniques applicable to satellite orbit determination, it is 
believed that they do not represent reportable or patentable items within .the 
meaning of the New Technology Clause.- Our reviews and final survey found no 
other items which could be considered reportable items under the New Technology 
Clause. 
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VII.' CONCLUSIONS AND RECOMMENDATIONS 


Simulations demonstrate that the J-Adaptive estimator with estimated 
noise statistics can automatically, without human intervention, and with a very 
simple orbit model (e.g. two-body model), estimate satellite orbits to an 
accuracy comparable with the data noise levels, when excellent, continuous 
tracking coverage is available. Such tracking coverage may be available from 
satelllte-to-satellite tracking. The fact that a two-body model can be utilized 
in the estimator prediction equation eliminates the need for numerical integra- 
tion of the orbit dynamics equations, which in turn results in a very fast 
(computational) orbit estimation algorithm. 

This estimator also tracks (estimates) the unmodeled acceleration vector. 
From these latter estimates a parametric model of the unmodeled accelerations 
may be recovered. The precision of tracking of these unmodeled accelerations 
is a function of the measurement noise levels (measurement precision) . Every- 
thing else being equal, the J-Adaptive estimator with a' priori "engineered” 
unmodeled acceleration statistics can track the unmodeled accelerations with 
higher precision than the estimator with estimated noise statistics. However, 
to engineer these statistics requires trial runs and simulations. 

In the absence of excellent tracking station coverage (single station pass), 
the J-Adaptive estimator is unsuccessful; there is a relatively long transient 
before the estimator locks onto the unmodeled acceleration. A consider version 
of the estimator ("adaptive consider mode") shows some promise in this situation, 
but further research is required in this area. 

On the basis of study results, the following recommendations are made: 

(1) The J-Adaptive estimator with estimated noise statistics should be 

applied to the satellite-to-satellite tracking problem. Simulations 
should be performed in this environment to verify estimator per- 
formance. This is a very natural application because of the great 
abundance of data available. This abundance of data is precisely 
where the J-Adaptive estimator is extremely efficient as opposed to 
batch estimation methods. 
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The J-Adaptive estimator should be studied in an environment con- 
taining more realistic measurement error models, including station 
location errors, timing errors, and so on. 

The single station pass tracking situation should be investigated, 
with a view, for example, to restricting uncertainties to selected 
subspaces of the state space. 

Other measurement types (angles) should be incorporated in J-Adaptive 
estimator studies. These may enhance the noise statistics estima- 
tion, thus improving tracking performance. 

Finally, parametric model recovery from the unmodeled acceleration 
estimates generated by the , J-Adaptive estimator should be studied.. 
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Appendix A - Oblateness Model ; 

A low order earth oblateness model utilized in the simulation studies is 
given below. In the equations which follow, R = is the satellite 

earth-fixed position vector expressed in km; r = [r|; y is the gravitational 
parameter of the earth; J^, and are zonal harmonic coefficients; C ^ , 

£> 22 ’ C^, ^31* ^33’ ^33 are tessera l harmonic coefficients; and 
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the acceleration due to oblateness, in earth-fixed coordinates and expressed 
in km/sec“ is given by 


AR 


OE 


= F + F + F 
J2 J3 J4 


- F - F 
22 


31 


- F 


33 


Numerical values for the coefficients above are 


u 

- 

3.986032 x -10 5 


(Ion 3 /sec 

J 2 

= 

1.08265 x 10“ 3 

4 

(km 2 ) 

J 3 

= 

-2.546 x 10“ 6 

4 

(km°) 

J 4 

= 

-1.649 x 10 -6 

4 

(km 4 ) 

'22 

- 

0.1536 x 10~ 5 

4 

(km 2 ) 

22 

= 

-0.872 x 10~ 6 

4 

(km 2 ) 

31 

=2 

0.2091 x 10" 5 

4 

(km 3 ) 

31 


0. 287 x 10" 6 

4 

(km 3 ) 

33 

S3 

0.782 x 10“ 7 

4 

(km 3 ) 

33 


0.226 x 10“ 6 , 

4 

(km 3 ) 


= 

6378.1641 

(km) 
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